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Abstract 

The motion of tiny heavy particles transported in a co-rotating vortex pair, with or without par- 
ticle inertia and sedimentation, is investigated. The dynamics of non-inertial sedimenting particles 
is shown to be chaotic, under the combined effect of gravity and of the circular displacement of the 
vortices. This phenomenon is very sensitive to particle inertia, if any. By using nearly hamiltonian 
. dynamical system theory for the particle motion equation written in the rotating reference frame, 

I one can show that small inertia terms of the particle motion equation strongly modify the Melnikov 

function of the homoclinic trajectories and heteroclinic cycles of the unperturbed system, as soon 
as the particle response time is of the order of the settling time (Froude number of order unity). 
The critical Froude number above which chaotic motion vanishes and a regular centrifugation takes 
place is obtained from this Melnikov analysis and compared to numerical simulations. Particles 
with a finite inertia, and in the absence of gravity, are not necessarily centrifugated away from 
the vortex system. Indeed, these particles can have various equilibrium positions in the rotating 
^Ljl reference frame, like the Lagrange points of celestial mechanics, according to whether their Stokes 

' number is smaller or larger than some critical value. An analytical stability analysis reveals that 

^ . two of these points are stable attracting points, so that permanent trapping can occur for inertial 

! particles injected in an isolated co-rotating vortex pair. Particle trapping is observed to persist 

c/3 ' when viscosity, and therefore vortex coalescence, is taken into account. Numerical experiments at 

large but finite Reynolds number show that particles can indeed be trapped temporarilly during 
vortex roll-up, and are eventually centrifugated away once vortex coalescence occurs. 
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1 Introduction - Particle motion equation 

CN , The motion of tiny particles in fluid flows has many unexpected features which have been studied 

' for decades in various contexts (chemical engineering, atmospheric dust, plankton transport in the 

^T) • ocean, planetesimal formation, etc.). Even the simplest particle transport model (that is passive 

, heavy non-interacting Stokes particles) coupled to any simple flow model (prescribed laminar flow) 

I lead, in general, to non-integrable differential equations with six degrees of freedom for the particle 

■ position and velocity components. The dynamics is therefore very rich and it is not surprising to 

observe complex motion emerge in a wide variety of natural flows where dropplets, solid grains 
or even biological particles are present. The non-integrability of the motion equation of these low- 
^ ■ Reynolds number particles is due to the gradients of the base flow, and this flow does not need to be 

. very complex for chaotic particle trajectories to emerge. Indeed, many simple fluid flows have been 

shown to transport tiny solid particles in a complex manner, and this motivated numerous theoretical 
or numerical studies (Maxey & Corrsin [14] ; Wang, Maxey, Burton & Stock [20] ; Mac Laughlin [11] 
; Fung [7j ; Tsega, Michaelides & Eschenazi [18j ; Rubin, Jones & Maxey [16j ; Druzhinin [5j ; Haller 
&; Sapsis [S]). Also, the understanding of the interaction between particles and elementary vortical 
structures could even help understand turbulent particle transport or turbulence modiflcation (see 
for example, in this spirit, Marcu et al. |I12]. Davila & Hunt [3], Ferrante & Elghobashi [6j). 

This paper deals with the motion of low Reynolds number heavy particles transported in a 
co-rotating vortex pair, i.e. the two-dimensional flow induced by two identical point vortices (flgure 
[1]). This is one of the simplest unsteady flows with a natural periodicity due to the mutual influence 
of the vortices. Gravity, acting in the (x, y) plane, is also taken into account. This choice has been 
motivated by the fact that various analyses of heavy particle motion in horizontal mixing layers. 
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submitted to subharmonic forcing, have shown that vortex pairing can significantly influence the 
structure of the particle cloud (Raju & Meiburg |15| . Chein & Chung [5], Kiger & Lasheras [10])- In 
particular, vortex pairing has been observed to increase particle dispersion and to enhance particle 
homogenization. Even though the flows investigated by these authors are much more complex than 
the one investigated here, we believe that the inviscid model used here can highlight and quantify 
some complex features of dust motion during vortex pairing. 

Firstly, we wish to show that non-inertial heavy particles injected in this flow can have chaotic 
trajectories under the combined effect of gravity and of the rotation of the vortices (section [2|) . This 
effect tends to increase particle mixing. Also, the effect of weak particle inertia on this phenomenon 
will be discussed (section [3|). Secondly, we will study inertial non-sedimenting particles, and show 
that some of them can be trapped by two attracting points rotating with the vortices (sectiond]). The 
former effect (chaotic particle motion) induces strong particle mixing, whereas the latter (particle 
accumulation) induces preferential particle concentration and could lead to aggregate formation. 
For the sake of simplicity, we will neglect the flow modiflcation due to the particles, as well as 
particle interactions. 

In this elementary flow one can easily show that the distance 2d between the vortices remains 
constant, and that the vortices rotate around the center point (here the point x = y = 0) with a 
constant angular velocity Q = T /Aird? , where F is the circulation of each vortex. Even though the 
flow Reynolds number T /u is large, the particle Reynolds number is assumed to be small, due to 
the small size of the inclusions. The simplest motion equation of such low- Reynolds number heavy 
particles, in the frame rotating with the vortices, is 

mp^^-^ = mpg + GiTfia^Vf - ^^^'^ - rup^e- rup^c (1) 

where a, nip and /x denote the particle radius, the particle mass and the fluid viscosity respectively. 
The term 7e = —Q.'^Xp is the acceleration of the rotating frame relative to the laboratory frame. The 
last term appearing in the motion equation is the inertial Coriolis force, where 7c = 20, x dXp/dt is 
the Coriolis acceleration. Gravity acts in the plane perpendicular to vorticity as indicated in Fig. 
[TJ In non-dimensional form, using d to normalize length scales and 1/17 to normalize time scales, 
we get (without renaming the variables) : 

where Xp{t) denotes the particle position at time t, t = Qmp/^Girna) is the non-dimensional 
response-time of the particle (Stokes number), and e*^ is the unit vector along the z axis. The 
vector ey{){t) = svcitex + coste^ is the upward vertical unit vector of the non- rotating frame. Vt 
is the non-dimensional terminal velocity and is assumed to be small throughout this paper (weakly 
sedimenting particles) : 

Vt = Tg/n^d < 1 

Also, it will be convenient to introduce the Froude number Fr = dQ'^/g = t/Vt- Because particles 
are much heavier than the fluid, Eq. ([2]) is valid even if the Stokes number r is not small, as added 
mass force, pressure gradient of the undisturbed flow, lift and Basset force are negligible. 

The fluid velocity fleld Vy, in the rotating frame, is steady and reads : 

/dip ^^p\ 

where the streamfunction is given by ijj(x,y) = — 21n|z^ — 1| -|- (x^ -|- y^)/2 (and z = x + iy) for 
the inviscid vortex pair considered here. The former term in ip is the flow induced by the two 
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vortices, and the latter corresponds to the opposite of the velocity of the rotating frame relative 
to the laboratory frame. The curl of Vj is equal to —26*2 everywhere, except at the vortex centres 
where it is infinite. The corresponding streamlines are plotted in Fig. [2]J^a) : they take the classical 
form of an eight shape around the location of the vortices (±1,0). This phase portrait has three 
hyperbolic saddle points located at (0,0) and A — (\/5, 0) and B = (— -v/S, 0). The point (0,0) 
has two homoclinic trajectories, whereas the two other points are related by a set of heteroclinic 
trajectories forming two heteroclinic cycles. Because this flow is two-dimensional and steady, no 
chaotic fluid point trajectories are expected to emerge. However, it will be shown in the next 
sections that the motion of particles has complex and unexpected features in such a flow. 

Three situations are examined. In section [2] we examine the case r = and < ^ 1, where 
inertia is negligible and a chaotic motion takes place under the effect of the unsteady term due to 
gravity in Eq. ([2]). Then we focus on the case < r ~ Vr 1 and examine the effect of weak 
particle inertia on chaotic sedimentation (section [3]). The case Vr = and r > is treated in 
section H] : it will be shown that non-sedimenting inertial particles have equilibrium positions in the 
rotating frame, two of which are stable. A numerical experiment, taking into account the viscous 
diffusion of the vortices, is shown in the last paragraph of section HI 



2 Chaotic sedimentation of inertia-free particles 

The case of non-inertial sedimenting particles can be addressed by assuming r — )• in the particle 
motion equation ([2]), so that all inertia forces vanish to leading order. In addition, we assume that 
Vt is small but non-zero. We then recover the simplest motion equation for heavy particles in fluid : 

ti X 

= Vf{X,) - VrSyoit) (3) 

which has been widely used so far (see for example Stommel |17|). The particle dynamics therefore 
corresponds to a steady two-dimensional flow (Vf ), plus a time-periodic perturbation due to gravity. 
The dynamical system ([3]), in the absence of perturbation [Vt = 0), has a homoclinic trajectory 
and a heteroclinic cycle : it is therefore tempting to check whether or not a finite perturbation 
(0 < *C 1) induces a homoclinic bifurcation, leading to chaotic aerosol motion in the vortex 
system. This point can be readily addressed by calculating the Melnikov function of any of these 
separatrices (say, Sj) : 

/+00 j<r 
X eyo{t + to)dt 

where to is the starting time of the Poincare section of the perturbed system (i.e. T-stroboscopy of 
the system with T = 2tt/Q) and Xo(t) is the solution of the undisturbed system corresponding to 
a point moving along the separatrix Sj (see for example Guckenheimer & Holmes [8J), that is : 

^ = VfiMt)) (4) 

If the Melnikov function has simple zeros when to varies, then the Poincare application of the system 
(or one of its iterates) corresponds to a horse-shoe map, and the particle dynamics is chaotic in the 
vicinity of Ej (Smale-Birklioff theorem). For the homoclinic trajectory Eg the coordinates XQ{t) and 
yo{t) of Xq are even and odd functions of time respectively (choosing Xq{0) equal to the intersection 
between Sq and Ox). We obtain, after removing the integral of odd functions : 



where 



Mo (to) = ao sin (to) 
«o = / \Xo{t) sint + yo{t) costjdt 

J —oo ^ ' 
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and xq and ijQ are the coordinates of clXo/dt. For the heterochnic trajectories Si and S2 the 
coordinates XQ{t) and yo{t) of Xq are odd and even functions of time respectively (choosing Xq{0) 
equal to the intersection between the separatrix and Oy.) We obtain, once again removing the 
integral of odd functions : 

Mi (to) = ai cos (to) 

where 




— xo{t) cost + yo(0 sin tjdt. 



for 1 = 1 and i = 2. In order to obtain the coordinates of Xo(t) we need to solve with 
the appropriate initial conditions specified above, on Sj. The equation of the separatrix Sj reads 
— 21n|2;^ — 1| + = ^(Sj), where r = \z\. By injecting this last equation into the undisturbed 
system ^ we are led to a purely numerical differential system which can be solved numerically to 
obtain Xo(t) and yo{t). From these the above integrals can be calculated and we get ao ~ —1.6, 
«! ~ 1.7 and 02 ~ —5.8. Hence, the Melnikov functions of the three separatrices have simple 
zeros. We therefore conclude that, as soon as gravity is taken into account, a bifurcation occurs 
in the particle dynamics and chaotic trajectories exist in the vicinity of the homoclinic trajectory 
So and of the heterochnic cycle Si U S2. Numerical solutions of equation jS]) confirm this result : 
Fig. [2] shows the trajectory of two particles injected respectively near the separatrix So and near 
the right-hand-side vortex. When Vr = these particles follow the streamlines, since the flow is 
steady. When Vr increases the particle injected near the separatrix has a very complex trajectory : 
it rotates in a chaotic manner around both vortices alternatively, as a consequence of the breaking 
of So. In contrast, the other particle, injected in a non-chaotic zone, has a regular trajectory around 
the same vortex. For higher Vy's the former particle wanders in several places of the domain. The 
Poincare sections of 20 particle trajectories, computed over 100 periods, are shown in figure [3l for 
Vt = 0.03. A stochastic zone is clearly visible in the vicinity of Sq and Si U S2, in agreement with 
the Melnikov theory presented in this paragraph. 



3 Effect of inertia on chaotic sedimentation 

When particle inertia, though small, is no longer negligible compared to gravity effects, that is if : 

< r ~ < 1 

i.e. the Froude number Fr = dVt^ /g is of order unity, then the particle motion equation ([2|) can be 
solved approximately by looking for an asymptotic solution of the form (Michael [4], Maxey ||13j ) 
dXp/dt = Vf{Xp) + tV^ + O(r^), and keeping Fr as a fixed parameter as r — )• 0. We are led to : 



dX, 



p 



dt 



1 



VfiXp) + r -- Cyoit) - VVf.Vf + Xp- 2e, x Vf\ + O(r^) (5) 



Like in the previous section we recover a non-autonomous dynamical system with two degrees of 
freedom (the phase space being the physical space), in the form of a perturbed hamiltonian system. 
The perturbation is now dissipative and contains the gravity term, the inertia term, the centrifugal 
and the Coriolis forces. The Melnikov function of a separatrix Sj now reads : 



Mi{to) = -— —jf X iyoit + to) dt + rrii 

^1' 7-00 



where 



I ^ dXo / d'^Xo , V? dXQ\ 
m.= I ^ ^x(-^+^o(t)-2e.x-^) (6) 
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and Xo{t), running on Sj, has exactly the same expression as in the previous section (solution of 
equation (H])). Hence, the calculation of the to-dependent part of the Melnikov function is straight- 
forward. For the separatrix Sq the Melnikov function reads : 

Mo{to) = sin (to) + m-o 

where ao < has been calculated in the previous section (for Si and S2 simply replace uq by at 
{i = 1, 2), mo by TTij (z = 1, 2) and sin by cos). The parameter rrii is a constant, and manifests the 
effect of particle inertia on the dynamics. The Melnikov function Mj(to) has simple zeros if 

Fr < ^ = Fri, (7) 

\'mi\ 

leading to chaotic particle dynamics in the vicinity of Sj. By injecting the numerical Xo(t) calculated 
above into Eq. ([B]), we get, for each separatrix : mo ~ —42, mi w —26 and m2 ~ 8. Because Fr, ao 
and ai are of order unity, we conclude that inertia has a huge effect on So and Si, as it tends to 
prevent the occurrence of simple zeros in their Melnikov function. Instead of the chaotic trajectories 
observed in the inertia- free limit, a regular centrifugation will take place in the vicinity of So and 
Si. The effect of inertia is less significant for S2. These calculations show that as soon as 

Fr > maxFrj = Fr2 ~ 0.7, 

i 

the stochastic layer around S2 will vanish also. 

In order to check this result we have solved numerically Eq. ([5]), for 1000 particles released 
slightly above S2 (near point C of Fig. EJa)), with r = 0.005 and Fr E [0.1,1.5] (Fig. i]). If 
a stochastic layer exists around S2, then some particles are likely to penetrate inside the zone 
bounded by S2 U Si. These particles are those who are located in the lobes formed between every 
two intersection point of the unstable manifold W^{B) (coming out of the hyperbolic point located 
near point B of Fig. [2)^a)) and the stable manifold Ty*(j4) (converging to the hyperbolic point located 
near A) . We observe on Fig. 3] that some particles indeed penetrate inside the zone bounded by 
S2US1 for Fr < Fr2. The detailed shape of the curve for Fr < Fr2 depends on the shape and position 
of the initial cloud. As soon as Fr > Fr2 however, the curve is flat and no particle penetrate into 
the zone bounded by Si U S2, in agreement with the Melnikov analysis. Indeed, when Fr > Fr2 the 
Melnikov function M2(to) is constant and non-zero, so that M^'^(A) and iy"(i?) never intersect, and 
particles located outside cannot penetrate into the zone bounded by S2 U Si. 



4 Trapping of non-sedimenting inertial particles 

When the particle response time is of order unity and gravity effects are small (Vy ^ t = 0(1)), 
particles are expected to be centrifugated away from the vortices. Here we show that this is not 
always the case, and that some particles, in spite of their finite inertia, can be trapped and remain 
in the vicinity of the vortex system. For the sake of simplicity we assume that the gravity term in 
the particle motion equation ([2]) can be neglected, so that the fluid velocity fleld in the rotating 
frame is now steady. By using this simplification we notice that Q has equilibrium points where 
the centrifugal pseudo-force balances the Stokes drag : 

^=Vf{Xp) + TXp. (8) 

If one of these points were asymptotically stable, which is the case if all the corresponding eigenvalues 
have strictly negative real parts, some particles could be trapped there and remain fixed in the 
rotating frame. This would imply permanent trapping by the vortex pair. In order to clarify this 
point we analyze, in the following lines, the equilibrium points and their stability. 
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The equilibrium equation has a trivial solution, Xp^^ = (0,0), for all r. By re-writing equation 
([8|) in polar coordinates {r,6), and after a few algebra (see appendix [A]), we obtain four more 
equilibrium points if < r < 2 — ^/S or r > 2 + ^/S. These points are symmetric with respect to 
(0,0) and are denoted ibXeg and ±Xeq . They correspond to 

4t 

sin 20 



r2 + l 



that is, 02 = arcsin(4r/(r2 + 1)) /2 (mod vr) and 6^ = tt/2 — 62 (mod vr). The coordinate r of these 
points is given by = (sin 29) jr + cos 20, that is : 



„2 



r 



r2 + l 



(4 ± v^r4 - 14r2 + 1) (9) 



The stability of these points can be addressed by re- writing the motion equation ([2]), with = 0, 
as a standard dynamical system with four degrees of freedom, with state variable Y = {x,y,x,y), 
where (x, y) denote the coordinates of the particle. Then ([2|) is an autonomous system of the 
form dY /dt = F{Y). One can check, by studying the roots of the characteristic polynomial of the 
Jacobian matrix VF = (dFi/dYj) at the equilibrium points, that ^Xeq is always unstable (see 
appendix [B]) . 

Let us now consider the equilibrium point Xgq . The determinant of VF at Xeq is : 

A(xg)) = -^(r'- Ur' + 1 - 4(r^ - Ur' + 1)'/') (10) 

It is negative only if — 14t2 -|- 1 > 16, that is r > VTS. Hence, by applying the same arguments as 
before, we conclude that for r > vTS this equilibrium point is also unstable. However, for smaller 
r, and in particular for all r in the range < r < 2 — ^/3, we have A{xiq^) > : the equilibrium 
point Xeq is not necessarily unstable. To prove that this point is indeed stable we need to solve for 
the characteristic polynomial of VF at Xgq , which reads : 

P(A) = a4 + ^ + (2 + 1)A2 + ^ + A 

T T'^ T 

where A(r) simply denotes A(Xeq ). By setting z = \ + l/(2r) we are led to a simpler polynomial : 
_ p 16r^z^ + 8t2(4t^ - l)z2 + 8T^(2r2A(r) - 1) + 1 

which can be readily solved. We obtain : 



2 _ l-4r2±4T2y^l- A(t) 

Because we consider the case A > 0, two cases emerge : < A < 1 and A > 1. In the former 
case (0 < A < 1) then z^. are real. Clearly, we have < 1/(4t2). If < then z_ is pure 
imaginary and the two corresponding A's have a negative real part equal to — 1/2t. If > then 
Z- < l/(2r), so A < 0. The second pair of roots are given by : 



4 - < V(4.^) 



since y^l — A(r) — 1 < 0. Here also we conclude that the corresponding eigenvalues of VF are real 
negative (or are complex with negative real parts if < 0). 
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We now turn to the case A > 1. By setting a = Re{z) we are led to o? = {Re{z^) + |z^|)/2 and : 

At'^o^ = (1 - + (1 - + 16r4A(r))i/2)/2 (11) 

One can prove that the function of r appearing at the rhs of this equation is always smaller than 
1 (see appendix [C]). Therefore, a < 1/(2t), and the corresponding eigenvalues of VF satisfy 
Re{X) = a - l/(2r) < 0. 

We therefore conclude that, for all r in the range < T < 2 - there exists a pair of 
asymptotically stable points located at zizXeq . 

In order to check these results, and to illustrate the effect of the attracting points ±Xeq on 
the particle dynamics, we have computed the trajectories of 10000 particles injected at t = 
in the square [—3,3]^, for r = 0.2. Results are shown in figure [5] where the particle cloud is 
plotted. Graph (a) shows that trapping indeed occurs, and that some particles converge to the 
points ibXeg predicted by the theory. Also, graph (b) shows the initial position of all particles such 
that \Xp lb Xeq \ < 0.2 at t = 20 : it gives an approximate view of a 2D cut, in the {x,y) plane, of 

the basin of attraction of the equilibrium points Xeq and —Xeq . 

In addition, we have conducted a series of runs, where 2000 particles were injected at t = 

all over the square [—3,3]^, for various response times r. We then have computed the percentage 

of particles lying inside the disk of radius 2 and centre (0,0) at long times (here t = 50). Most of 

these particles are those which have been "trapped" by the two attracting points, if any. Clearly, 

in the absence of attracting points p{t) ~ 0, whereas p(t) is finite if attracting points are present. 

So, if T is not too small this percentage is a good indicator of the occurrence of trapping. Figure [6] 

shows that, as soon as r > 2 — \/3, the percentage is identically 0, in agreement with the stability 

^(3) 

analysis : the determinant is negative there, so that the point iXeq is unstable. 

Comparison with a numerical experiment. These results are valid for inviscid fluids. When viscosity 
is finite the flow is significantly different since vortices coalesce and finally form a single vortex. 
Particles will therefore be centrifugated away once coalescence occurs. However, if the diffusive 
time scale over d is larger than the turnover time one can expect particles to be trapped 

temporarilly in the vicinity of ^Xeq (the stable points of the inviscid flow investigated above) during 
the interaction of the vortices, even though the equilibrium points do not exist strictly speaking 
since the relative flow is no longer steady. We have therefore performed a series of runs with a 
two-dimensional (spectral) Navier-Stokes solver including a lagrangian particle tracking algorithm, 
within a periodic box [—ird, nd]'^. The initial condition corresponds to a pair of Lamb-Oseen vortices 
Vf{x,y) = {-u{r)y,uj{r){x - x^)) with 

w = r[l - exp(-3rV'^^)]/(2vrr2) 

where = (x — x^)^ + y^ and = ibd is the initial x-coordinate of the vortices. The circulation of a 
single vortex of this kind, placed in an infinite domain, is independent of the core thickness 5 and is 
equal to F. The fluid and particle motion equations are set non-dimensional by using = AirtP /T 
for times and d for scales, like for the inviscid calculations presented above. The initial core-size 5 of 
the vortices is taken to be one tenth of the flow domain. The Reynolds number Re = Qtf' /v is equal 
to 800. The number of Fourier modes is 256 in each direction, and a second-order Adams-Bashforth 
algorithm is used for time stepping. Particles are passive (they do not modify the flow), and do 
not interact. They are injected at random at i = and cover the whole vortex pair. Their initial 
velocity is equal to that of the fluid, and their trajectories are calculated by solving equation Q 
(without the inertia forces since the numerical solver used here considers the dynamical equations 
in the laboratory reference frame). 
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Figure [7] shows the evolution of the vortices and of the particle cloud, when r = 0.1 (this value 
has been chosen because it corresponds to the peak of figure [H]). We observe that particle trapping 
indeed occurs, and persists until t ~ 12 (that is about two turns). Because this calculation is done 
in a periodic domain, the influence of the neighbouring vortex pair located in boxes [— tt + 2n7r, tt + 
2n7r]^, n S Z is signiflcant here. However, these vortices do not prevent the appearance of the 
attracting points predicted by the inviscid theory in inflnite domain (section . When coalescence 
begins, the set of trapped particles is elongated into a thin filament (Fig. [7|) which will then be 
centrifugated away. The dashed line in figure [7] is the trajectory of the relative equilibrium points 
ibXeg of the inviscid theory (equation ([9])). This line is close to the position of the accumulation 
points, in spite of the effect of viscosity and of the neighbouring vortices. 

5 Conclusion and open questions 

We have investigated the motion of tiny heavy particles with or without inertia and gravity in a co- 
rotating vortex pair. For sedimenting inertia-free particles (r = and < ^ 1) we have shown 
that, under the combined effect of gravity and of the rotation of the vortices, a chaotic particle 
dynamics can take place. By writing the particle motion equation in the rotating reference frame 
attached to the vortices, we observed that the particle dynamics take the form of a hamiltonian 
system submitted to a periodic forcing due to gravity. A "stretch-sediment and fold" mechanism is 
responsible for this chaotic particle dynamics, like the one described in Ref. [Ij, as gravity plays a 
central role in the folding of volume elements of the dispersed phase. 

Chaotic motion is very sensitive to particle inertia when < r ~ ^ 1- Indeed, by using 
nearly hamiltonian dynamical system theory for the particle motion equation written in the rotating 
reference frame, we have shown that small inertia terms of the particle motion equation strongly 
modify the Melnikov function of the homoclinic trajectory of the unperturbed system. In particular, 
the stochastic layer in the vicinity of the separatrices Sj (i=0,l,2) vanishes as soon as the Froude 
number is larger than some critical values Fr^ given by equation jT]). A regular centrifugation 
therefore takes place as soon as the Froude number is above maxFr^ ~ 0.7. Numerical results 
confirm this value (figure U]) . 

Particles with a finite inertia, and in the absence of gravity (r = 0(1) and Vr = 0), can have 
various equilibrium positions in the rotating reference frame, according to whether their Stokes 
number is smaller or larger than some critical values of order unity. We have rigorously shown 
that two of these points are stable attracting points, so that permanent trapping occurs for inertial 
particles injected in an isolated co-rotating vortex pair. Numerical computations confirm this result, 
and show that the basin of attraction of the attracting points covers a non-negligible part of the 
flow domain. Our analytical calculations also show that trapping should stop if r exits the range 
< T < 2 — \/3. This point has been confirmed by computing the percentage of trapped particles 
(which can be thought of as a measure of the basin of attraction) , for various reponse times r : this 
percentage drops to as soon as r > 2 — -y/S. Note that the points ^xiq'^ and ^xiq^ also exist for 
r > 2 + "v/S- However, the former has been shown to be unstable for all r, and the latter is also 
unstable as soon as r > \/T5- Numerical experiments taking into account the effect of viscosity 
have been conducted. We observe that, when the Reynolds number Qcf /v is large, the effects of 
viscosity are sufficiently slow to enable particle trapping in the vicinity of two points, the position 
of which is in qualitative agreement with the results of the inviscid theory. Once vortex coalescence 
is complete, particles are centrifugated away. 

Note that this effect could be of interest also in the context of particle motion in protosellar 
gaseous disks. Indeed, large scale anticyclonic vortices might play a key role in dust trapping, with 
the help of the Coriolis force due to the disk rotation, and these vortices are known to interact and 
coalesce. The dynamics of dust during vortex pairing in protostellar disks could therefore be a topic 
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of interest. 

When r = 0(1) and Vp is non-zero the flow equation is unsteady in the rotating frame, and the 
attracting equihbrium points described above no longer exist. The analytical treatment of particle 
dynamics in this case is much more complex, but numerical simulations, not shown here, suggest 
that trapping could exist anyway. Further analyses should clarify this point. 

The generalization of these results to vortex pairs with different strengths is the next step of this 
work. Preliminary numerical simulations show that particle trapping still exists in the asymetric 
case, but the analytical treatment is heavier and is currently under study. Also, when vortices are 
located in the vicinity of a wall and move parallel to it, particle trapping is expected to persist. 
Indeed, Vilela & Motter |19| have observed attracting points in their simulation of particle transport 
in a leapfrogging vortex system (which also corresponds to a 2D inviscid vortex pair in the vicinity 
of a wall). We therefore conclude that the attracting points persist in the presence of the wall, 
even though these points have a more complicated trajectory. A detailed theoretical analysis would 
enable to determine the range of parameters leading to such a trapping. 
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A Equilibrium positions of inertial particles 



The equilibrium equation ([S]) reads 

-4y(r2 + 1) 



r4 + 1 + 2(y2 _ 2;2) 
4x{r^ - 1) 



y + TX = 0, (12) 
x + ry = (13) 



+ 1 + 2(y2 _ 2;2) 

where = + gy using polar coordinates, and combining these two equations, we are led to : 

= (sin26l)/r + COS26', (14) 
together with the admissibility condition required to avoid a null denominator : 

+ 1 - 2r2cos26' / (15) 
By injecting this r into ()12p -(|13 p we get : 

sin 29 (t sin 9 — cos 9) [(l + r2)sin26i-4r] =0 

and 

sin 29 (sin 9 + t cos 9) [(1 + r^) sin 29 - 4r] =0 

Because r sin 9 — cos 9 and sin 9-i-T cos 9 cannot be both zero, two families of solutions emerge. The 
first one is sin 20 = and r = 1 : it does not satisfy the admissibility condition (jlSp . and must 
be rejected. The second one exists if 1 + > 4r, that is < r < 2 — ^/S or r > 2 + ^/S, and 
corresponds to : 

4t 

sin 261 



1 + t2' 

that is ^0 = arcsin( ^^^^ )/2, and 9 = it + 9q, and 9 = 7r/2 — 9o, and 9 = 37r/2 — 9o (mod vr). These 

^(2) ""(3) 

angles corresponds to 4 equilibrium positions in the physical plane, denoted ±Xeq and ±Xeq . The 
r coordinates of these points are obtained by using ([1] 



B Instability of equilibrium point (2). 

The determinant A{xif) of VF at X^g'' is : 

A{xg>) = (r^ - 14r2 + 1 + 4{t' - Ur' + 1)'/') 

and one can check, since 2 it -v/3 are also roots of — 14r2 + 1, that — 14r2 + 1 is strictly positive 
for all r such that < r < 2 — ^/S or r > 2 + y^. We then conclude that A{xig^ ) is always negative 
in this range of r. This implies that the equilibrium point Xeq is unstable. 

->(2) 

Indeed, let Aj, (i = 1, 4) be the roots of the characteristic polynomial of VF at Xeq . Because 
the coefficients of the polynomial are real these roots are either real or complex conjugate. Also 
A1A2A3A4 = A{Xeq ) < 0. Clearly, this implies that the four roots cannot be complex conjugates 
(i.e. A2 = Ai and A4 = A3) since we would have A{xig^) = |Aip IA3P > there. If two of these roots 
(say Ai and A2) are complex conjugates and the other two are real, then we must have IA1I2A3A4 < 0, 
so A3A4 < : hence there exists a strictly positive eigenvalue : the equilibrium point is unstable. 
Finally, if all the roots are real, then, their product being strictly negative, one of them at least is 
positive : Xeq is unstable. 
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C Upper bound for Re{zY 

To prove that a? < l/(4r^) we need to show that (Eq. (jlip ) : 

1 - + (1 - 8t^ + 16T^A(r))^/2 ^ 2 

for ah r in the range < r < 2 — \/3- This is equivalent to : 

(1 - 8t2 + 16r^A(r))i/2 < 1 + 4^2 

Taking the square of these positive numbers, and after simplifications, we are led to : 

T^A < 1 + 

From (fTO|) we get A = (—q + 4y^)/(4r^), with g = — 14r^ + 1. The last inequality is therefore 
equivalent to : 

Once again, taking the square of these positive numbers, we are led to an equivalent expression : 

(7 - (| + 1 + r2)2 < 

Because q = — 14t^ + 1 we have: 

g - (| + 1 + r^)^ = [-9 - 124r2 - r^(94 - 20t^) - t^]/16 

and this last quantity, being a sum of negative numbers (since < r < 2 — VS), is negative. We 
therefore conclude that < l/(4r^). 
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Sketch of the two vortices (black dots), together with the rotating reference frame 
with axes (0,x,y) and the laboratory frame with axes (O,xo,yo)- The lower graph 



is the homoclinic trajectory of the fluid points dynamics in the rotating frame. ... [14 

2 Streamlines in the rotating frame (a). Trajectories of two inertia- free particles for 
various terminal velocities (b)(c)(d). The particle injected near the vortex has a 
regular trajectory around it. In contrast, the particle injected near the separatrix Sq 
wanders in various places of the flow domain 

3 Poincare section of 20 inertia-free particles, for Vr = 0.03 

4 Plot of the percentage of particles released at t = slightly above S2, and crossing S2 
during the simulation (the final time is t = 100). As soon as Fr > Fr2 the stochastic 
layer located in the vicinity of E2 vanishes and no particle crosses E2, in agreement 
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with the Melnikov analysis [17 

(a) : plot of a particle cloud after 20 convective times. At t = the cloud covers 
the square [—3,3]^, and is composed of inertial particles with r = 0.2 and = 0. 
The equilibrium points X^q are indicated by the white dots (i). Graph (b) shows the 
initial position of trapped particles [3 
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Plot of the percentage of "trapped" particles at t = 50, versus r, when Vr = 0. As 
soon as r > 2 — -v/3 ~ 0.27 no trapping is observed (p(r) =0), in agreement with the 

stability analysis 

Evolution of a particle cloud (in the laboratory frame) advected by a vortex pair 
at large but finite Reynolds number, for r = 0.1 and Vr = 0- Lines are vorticity 
contours. Upper figure : t = 10, lower figure : t = 15 
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Figure 1: Sketch of the two vortices (black dots), together with the rotating reference frame with 
axes (0,x,y) and the laboratory frame with axes {O,xo,yo). The lower graph is the homoclinic 
trajectory of the fluid points dynamics in the rotating frame. 
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-3-2-10123 -3-2-10123 

Figure 2: Streamlines in the rotating frame (a). Trajectories of two inertia- free particles for various 
terminal velocities (b)(c)(d). The particle injected near the vortex has a regular trajectory around 
it. In contrast, the particle injected near the separatrix Sq wanders in various places of the flow 
domain. 
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Figure 3: Poincare section of 20 inertia-free particles, for Vr = 0.03. 
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Figure 4: Plot of the percentage of particles released at t = slightly above S2, and crossing S2 
during the simulation (the final time is t = 100). As soon as Fr > Fr2 the stochastic layer located 
in the vicinity of S2 vanishes and no particle crosses S2, in agreement with the Melnikov analysis. 
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(a) y (b) 




-3-2-10123 -3-2-10123 



Figure 5: (a) : plot of a particle cloud after 20 convective times. At t = the cloud covers the 
square [—3,3]^, and is composed of inertial particles with r = 0.2 and Vr = 0. The equilibrium 
points Xeq are indicated by the white dots (i). Graph (b) shows the initial position of trapped 
particles. 
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Figure 6: Plot of the percentage of "trapped" particles at t = 50, versus r, when Vt = 0. As soon 
as r > 2 — "v/S ~ 0.27 no trapping is observed (p(t) = 0), in agreement with the stability analysis. 
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Figure 7: Evolution of a particle cloud (in the laboratory frame) advected by a vortex pair at large 
but finite Reynolds number, for r = 0.1 and Vr = 0. Lines are vorticity contours. Upper figure : 
t = 10, lower figure : t = 15. 
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